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Abstract 

We investigate whether the equilibrium time averaged state of a self-organising 
system with many internal degrees of freedom, 2D- daisyworld, can be described 
by optimising a single quantity. Unlike physical systems where a principle of maxi- 
mum energy production has been observed, Daisyworld follows evolutionary dynamics 
rather than Hamiltonian dynamics. We find that this is sufficient to invalidate the 
MEP principle, finding instead a different principle, that the system self-organises to 
a state which maximises the amount of life. 
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I. INTRODUCTION 



Thermodynamics provides an excellent macroscopic description of the physics of matter 
in terms of time-independent averaged quantities such as energy, entropy, free energy. These 
arise from the underlying microscopic, time-dependent motions of atoms. The two descrip- 
tions are unified through statistical mechanics, and therefore contain the same information, 
but in in the absence of atomic-level detail the thermodynamic picture is far more useful for 
everyday application. 

It is interesting to consider whether ecosystems can similarly be described by macroscopic 
averages rather than details of individual species. 

Recent work has suggested that the principle of maximum entropy production (MEP) 
may be a good way to investigate complex, self organising systems. Lorenz et al (2001,2002) 
suggest that the atmosphere of Earth and Titan may be in such a state. Dewar (2003) has 
shown that other known distributions of open, dynamical systems, such as self-organised 
criticality (Bak et al, 1988), the fluctuation theorem (Crooks, 1999) and Zipf's Law (Zipf 
1932). The basis of Dewar's argument is that if one takes an appropriate average over all 
its internal possibilities Q], consistent with the known fluxes and internal energies, MEP will 
emerge. Of course, real systems do not actually sample all possibilities, so this only provides 
a sufficient conditions for MEP, the necessary condition is that the average over possibilities 
actually observed converges rapidly. This parallels an argument in equilibrium statistical 
mechanics, where thermodynamic averages are calculated from an ergodic hypothesis (that 
all microstates contribute according to Boltzmann statistics) to which the chaotic trajectories 
of real systems rapidly converge. 

Boltzmann statistics follow from Hamiltonian dynamics. In the case of ecological systems, 
the dynamics are driven by evolution, and no proof of convergence of an evolving system 
to a Boltzmann-type distribution yet exists. It is interesting to ask whether more general 
principles, such as MEP, can be applied in an ecological context. In this paper we explore 
a particular ecological model, 2D-Daisyworld, in which the MEP state can be rigorously 
defined, to see whether the system does, in fact, self-organise thereto. 

The 2D-Daisyworld, is based on previous work (Von Bloh et al., 1997,1999, Ackland et 
al 2003). Life is reduced to a single species 'daisies' with a single property (albedo) upon 
which selection can act. The local environment is reduced to a single variable temperature. 
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The planet is warmed by a Sun whose heat production (insolation) increases slowly. Growth 
is only possible across a narrow range of temperatures (5-40''C) with the mid point 22.5°C 
being optimal and the growth rate dropping parabolically to zero at the extremes. Local 
temperature is partly determined by insolation and albedo and partly by diffusive heat flow 
across the planet. Growth occurs by seeding of bare ground adjacent to existing daisy 
populations, eliminating the need in the original, zero-dimensional daisyworld of Watson 
and Lovelock, (1983) to predefine different local environments for different albedo types. 
Each (asexual) offspring mutates in albedo at random. 

The original daisyworld is deterministic, cast in the form of differential equations, and has 
spawned many variants to incorporate specific evolutionary dynamics or natural selection 
(Lenton, 1998, Lenton and Lovelock 2000, Lovelock 1998, Robertson and Robinson, 1998, 
Staley 2002, Sugimoto, 2002). The 2D version is stochastic, with many internal degrees 
of freedom, allowing an entropy production to be defined by the range of albedo (called 
biodiversity), as well as through heat flow. 

Historically, daisyworlds were designed to illustrate the self-regulating effects of couphng 
between life and the environment, and their impact on evolution by group selection. This 
is not the issue of interest here - daisyworld is used as a test case for MEP. 

II. METHODS 

The temperature field changes as: 

Ct{x,y,t) = DTV''T{x,y,t) 

+SL{1-A{x,y,t)) (1) 

Where C = 2500 is the heat capacity, T the temperature field, Dt = D/C the thermal 
diffusion constant = 5.67 x 10^^ the Stcfan-Boltzmann constant, Sl = S'/432.3, the 
insolation normalised to that which gives T=:295.5K on a bare planet and A the albedo 
field (daisy/bare ground). Each lattice point has an associated temperature and albedo. A 
ranges from to 1, with 0.5 being the value for bare ground. The daisy field is discrete and 
evolves stochastically. We use a square lattice with eight neighbouring sites. At each time 
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step each site is examined and: 

(1) If occupied by a daisy, it changes to bare ground with probabihty 7(T) (the "death 
rate"). 

(2) If occupied by bare ground, the site is populated with probability /3(T) (the "growth 
rate"). 

The new daisy has the albedo of its parent, from a randomly chosen neighbouring site, 
with a random fluctuation (mutation) drawn from a uniform distribution between r and -r. 
If the neighbouring site is bare, no growth occurs. 

III. MEP OR GAIA? 

Two global principles will be advanced, these are not new, but formal definitions for the 
present system are made. 

1/ Maximal Entropy Production (Dewar, 2003) - the system self organises to produce 
entropy as rapidly as possible 

2/ The Gaia principle (Lovelock, 1998) - the system self organises to maximise the amount 
of life 

In previous work on the 2D case (Ackland et al 2003), it was shown that multiple solutions 
satisfied the Gaia principle. A secondary criterion was found, maximisation of biodiversity 
entropyP] A^^lnA^ where is the number of daisies with albedo A, leading to a mean 
distribution 

{Na) = e^^ (2) 

with (3 given by the solution of: 

{A) [\xp{(3A)dA = [\exp{(3A)dA (3) 
Jo Jo 

Thus the global state of 2D daisyworld was classified in terms of two quantities, maximised 
growth (i.e. temperature regulation at the favoured growth temperature) and maximised 
biodiversity. These are equivalent to satisfying both MEP and Gaia principles simultane- 
ously. 

It is possible to make a small change to the 2D daisyworld so that Gaia and MEP cannot 
be satisfied simultaneously. To do this simply requires making the death rate temperature 
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dependent, and the temperature for minimum death rate (T^) different from that for max- 
imum birthrate {Tg). In this case the MEP will predict that the temperature will adjust 
close to Tg, while the Gaia hypothesis predicts it will adjust close to T^. The explanation 
for this is straightforward: MEP predicts that entropy (biodiversity) be produced at the 
fastest possible rate: maximising the birthrate. By contrast, the Gaia hypothesis is that the 
total number of daisies be maximised: minimising deathrate. Thus investigating whether 
the mean temperature tends to Tg or Td determines which (if either) of the two hypotheses 
applies. Figure^ shows the variation of mean temperature and population as a function of 
insolation: the mean temperature obtained is and the population is 1 — 7(Td). Thus we 
conclude that the dynamics satisfy Gaia rather than MEP. 

IV. THERMAL ENTROPY PRODUCTION 

The arguments above pertain to biodiversity entropy. With a finite diffusivity, one can 
also investigate whether the MEP can be applied to entropy production by heatflow in 
daisyworld. There are conceptual difficulties with defining the system here - the models 
for MEP on Earth (Ozawa et al, 2003), and Dewar's (2002) formulation assume constant 
fiows of energy into and out of the system; the self-organising albedo of daisyworld is able 
to regulate these fiows. Similarly, the diffusion constant is fixed in 2D-daisyworld, whereas 
in atmospheric models it is assumed to self adjust. Notwithstanding this, there is a range 
of temperatures in 2D daisyworld, and heat is transported from hotter to colder regions, so 
we can define thermal entropy productionjj] as 



Given the constraints on fiow, it is not straightforward to determine the maximum possi- 
ble entropy production. However, after an abrupt change in the external forcing (insolation) 
the system will be out of equilibrium, and will then evolve towards a steady state. If MEP 
holds, the state immediately after the change should have lower entropy production than 
the steady state. Figure El shows no sign of such behaviour. Monitoring entropy production 
also enables us to monitor another property, the net heat fiow in and out of the system. In 
steady state, the mean values of these heat fiows should be equal - in OD daisyworlds solu- 
tions initially assumed input and output were in balance at all times - this is not required 




i k=l 



(4) 
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for a nonlinear system, but proof of the stability of these solutions came later (Saunders, 
1994). In the 2D-daisyworld Figure El shows the time variation of flow in and flow out. An 
interesting feature here is that flow in (basically albedo) fluctuates much faster than flow 
out (basically temperature), which also has a delayed response of many daisy lifetimes. 

Observations of planetary MEP are presumed to indicate that large-scale atmospheric 
heat transfer mechanisms are self-organising. However, even in the curved-planet version 
of 2D-daisyworld (Ackland et al, 2003), very strong temperature regulation eliminates large 
scale temperature variation and hence the driving force for the atmosphere (except in the 
case of desert formation). Figure El shows that the mean heat entropy produced per cell 
does not tend to a maximum during equilibration after a sharp change in insolation. We 
this conclude that no MEP principle applies to heatflow in 2D daisyworld. 

We examined whether there was an optimum value for the diffusion constant D: at given 
S, heat entropy is largely independent of D. This is non-trivial, since it implies that for high 
D temperature variations are suppressed precisely enough to compensate for the increased 
diffusion constant. However, it indicates that even if D could be optimised by a hypotheti- 
cal daisyworld atmosphere, in the presence of strong temperature regulation no detectable 
maximum entropy producing value of D exists. 

V. MEP IN ORIGINAL DAISYWORLD 

In Watson and Lovelock's (1983) OD daisyworld two daisy albedos exist, denoted by 
suffices b and w, with albedo Af, = 0.25 and A^j = 0.75 respectively. The relevant equations 
are: 

Change of fractional areas occupied by each daisy type (set to zero at equilibrium) 



dai / dt 



ai{x(3i - 7); 



(5) 



Fraction of free surface for growth 



X = 



1 - Zl^i! 



(6) 



Growth rates for daisies. 



A = max{0, 1 - ((22.5 - Ti)/17.5)2}; 



(7) 
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Mean planetary albedo (reflectivity) 



Ap — x/2 + Aitti] 



(8) 



Local temperatures over daisyfields 



(T, + 273)^ = q{Ap - A,) + (Te + 273)^ 



(9) 



One further equation closes the set - this represents the equality of absorbed and emitted 
heat, As we saw in the 2D case, this equality strictly should be considered as a time average, 
however if we assume that no variable large heat reservoirs (e.g. icecaps) exist the timescales 
are such that the following equality holds: 



S being an insolation, L a variable of order unity, and emitted by Stefan-Boltzmann 
radiation. 

The usual conclusion from this model is that the daisies moderate the temperature, 
however the graph of T f s L (figure 0^): shows that Tg is not optimised to the optimum for 
growth. In contrast, a plot of x vs L (figure |31d) shows that x is a minimum for all L. 

The daisyworld equations can be solved analytically. The solutions fall into four inde- 
pendent classes: no daisies, black-only, white-only, both types. At a given insolation only 
one "living" solution is stable if the simulation is started with non-zero concentration of 
both daisy types. Dead and living solutions coexist in regions where the dead planet tem- 
perature lies outside the limits for growth (5-40°C). The stable solutions (FigE)) have been 
found dynamically, (iterating equations IMTUl to self consistency) by numerous authors. The 
important point here is that the stable solution can be determined by minimisation of x 
among allowed solutions. Adopting the lowest x is the OD equivalent to the Gaian principle 
of maximum life, as found in the 2D case. 

As a further test, we introduce a third daisy species (grey) with albedo 0.5, the same 
as bare ground and therefore with no regulating effect. Equations El show that coexistence 
of three species is impossible, and in fact the stable solutions (Figj^J traverse regions of 
black-only, black and grey, grey only, grey and white, white only. In each region unstable 
equilibria exist for other combinations, and in each case the unstable equilibria have higher 




(10) 
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X than the total one. This result is robust against changes to g, Ai and 7 (including species 
dependent 7). 

The equilibrium solution to daisyworld can also be obtained by either integrating equation 
El or with a minimisation of x with respect to subject to constraints. Each gives the same 
solution f Figs mis]) , independent of the history of the system. Thus the equilibrium state 
in 2D and OD daisyworld can be determined assuming the Gaia principle, without knowing 
detailed dynamics, just as the equilibrium states of thermodynamic systems can be determined 
by minimising free energy. 

VI. GAIA AND MEP IN THE LOGISTIC MAP 

The logistic map is another widely used model for population growth (Hassell, Lawton 
and May, 1976). The resulting populations famously follow the period doubling route to 
chaos. There is no natural selection and again it is interesting to ask whether this model 
displays any kind of maximum population/maximum entropy behaviour. 

The classical logistic map gives an nth population size of 

y{n) = ay {n - - y in - I)] (11) 

The equilibrium distribution for y depends only on the parameter a. The mean popu- 
lation, averaged over 40000 iterations, is shown as a function of a in fi^ It exhibits a lot 
of structure which can be connected with the bifurcation diagram of the logistic map. The 
mean population of the attractor is higher than would be expected from maximum entropy 
(all levels equally occupied). 

Turning now to entropy, there are many definitions related to the population distribution 
by analogy with the 2D daisyworld. For example, the Tsallis entropy (Latora et al 2000) 

N 

S = J2 log{P{yi))/Nlog{M/N) (12) 

where P{yi) is the number of occurences of a population between ^ and during M 
iterations of the logistic map (taken after a 1000 iteration equilibration period). 

This quantity has a maximum possible value of 1, and (neglecting numerous distinct 
special cases of limit cycles) its value tends relatively smoothly to this value with increasing 
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a. However, its actual value depends on N. This only tells us that the area lying within the 
attr actor increases with a. 

The attractors in the logistic map provide some constraints on the mean value of y and 
the Tsallis entropy: it is difficult to define what the "maximum life" solution would be. 
Thus while it is notable that perturbed logistic map dynamics gives mean population in its 
stable region above 0.5 for all a, no quantitative principle emerges. 

VII. CONCLUSIONS 

The 2D-daisyworld model has been used as a testbed to determine whether any global 
maximisation is applicable to a system governed by the dynamics of natural selection. We 
find that Maximum Entropy Production is not valid here, either for biodiversity or heat 
transfer entropy. This casts doubt on some recent work in which the MEP is assumed as a 
constraint and used to obtain a solution (e.g. Pujol, 2003). Analysis of 2D and OD cases 
suggests that daisy worlds appear to obey a separate optimisation principle, in the simple 
2D case indistinguishable from MEP, that of maximising the amount of life present in the 
system. The logistic map attractors have above average population, but for this system it 
is difficult to define what the maximal population value should be. 

The more fundamental implication is that stability of a dynamic equilibrium can be 
determined as an optimisation problem on the properties of the system as observed, and that 
detailed integration of system dynamics through time may be unnecessary in determining 
the equilibrium state. 
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[1] This is implicit in using the sum over all paths to build a partition function 
[2] Finite difFusivity leads to clump formation, which confounds this result due to neglect the 
entropy associated with variation in clump size. In principle this could be incorporated, but for 
the present purposes the maximum entropy principle observed in the limit of infinite diffusivity 
is sufficient. 

[3] We say "close to" : these are the limiting cases for low death rates and high birth rates, assuming 

uncorrelated dead patches. 
[4] factor of two arises from double counting 
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FIG. 1: Fractional daisy coverage and temperature (in units of optimal growth temperature, 
Tg = 295.5 ) from a calculation with a 200x200 lattice (7 = 0.02(1 - {5 + Td-T){T + 6 - To) 126'^, 
P{T) = {S + Tg-T){T + S-Tg)/6^. = 305.5, 6 = 17.5)D = 0.1, r = 0.02) with scaled insolation 
S varied from 0.7 to 1.7 over 10^ updates of daisy and temperature fields with a temperature 
dependent death rate, minimal at = 305.5. The graph shows that in the non-desert regime 
population tends to its maximal value (1 — a; = 7/(7 + /3) = 0.99 in the mean field approximation, 
where regrowth is always possible) while temperature is regulated at T^, rather than Tg which 
would maximise(bio-)entropy production. Thus the governing global principal appears to be Gaian 
(maximal life) rather then MEP. At higher S, fluctuations lead to long lived deserts, which cause 
the mean field approximation to break down as daisies cannot grow in the interior. 
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Time (Mean daisy lifetimes) 



FIG. 2: Graph showing net total energy flows into (solid line) and out of (dotted line) a 200x200 cell 
2D daisyworld system (diffusion = 0.1, death = 0.02, mutation = 0.02). The simulation starts with 
S=1.0 and increases by discrete jumps of 0.1 every 10000 timesteps. Up to S=1.7 this is sufficient 
to allow the system to return to equilibrium: the excess incoming flow of energy is reduced by 
evolution of the albedo. Beyond 1.7 (see inset) the sudden change generates sufficient die back 
to create deserts in the system, which initially cause positive feedback and increase the incoming 
flow still further. Ultimately, evolution of paler daisies in the non-desert regions allow reinvasion 
of the deserts, but the simulation had to be run for 20000 timesteps (S=1.7) and 50000 (S=1.8) 
to achieve this. Even in equilibrium, transient desert regions occur. At S=1.9 the simulation was 
run for 120000 timesteps, throughout which desert and populated regions coexisted and migrated 
around the planet (there being no curvature to favour one region over another). Throughout, note 
that the energy emitted is smoother than that absorbed and lags the energy absorbed by about 
500 timesteps (10 mean daisy lifetimes) 7 = 0.02, D = 0.1,r = 0.02. 
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Time (Mean Daisy Lifetime) 

FIG. 3: Graphs showing variation of configuration (biodiversity) entropy {J2'i=iP{Ai) In p{Ai)) 
from partitioning the albedos into 20 groups, and (heat) entropy production (/ dQ/T) averaged 
over the 200x200 lattice and over a period of 50 timesteps (one mean daisy lifetime) for the same 
calculation as Fig|21 Each changes monotonically as the temperature is increased through the 
region of full coverage. Beyond S=1.7, where fluctuating deserts appear, each measure becomes 
less well defined. There is significantly increased heat flow from the hot deserts to the covered 
regions, so the heat flow becomes dominated by the area of desert. Conversely, the biodiversity is 
reduced when much of the planet is unpopulated. Reversing the change in temperature shows that 
the entropy change is reversible. 
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FIG. 4: Results from OD daisyworld with parameter values from Lenton and Lovelock: S/ctb = 
1.68 X l^y^^K^^ q = 2.06425 x 10^ and 7 = 0.3. a) Temperature variation with increasing solar 
insolation. Thick line show planetary average, dotted lines show Tj, and Ty^. The four regimes are 
shown in their region of stability: Simulations commencing with temperatures outside the growing 
range can stabilise the bare planet solution, but single-daisy and two-daisy solutions do not coexist. 
Only stable solutions are shown, including the bistable living/dead regime, b) coverage: Thick line 
denotes bare ground, thin lines denote individual species. For non-zero at and a^, da-i/dt = 
means that /3b = f5uj and from [7| and IHl Tyj and T;, are also independent of L: = 22.5 for q = 
falling almost linearly to zero for q = 9.34 x 10^. For standard parameter values from Watson 
and Lovelock, (3b = = 0.918 Tw=17.5 and Tf,=27.5. With and (3^ fixed by the choice of q, 
independent of L, as is the amount of bare ground at equilibrium: x = "f/ fiw] (from 1). Equations 
and El now give us simple linear expressions fey ab and in terms of x and the single variable 
Ap. Thus, once the assumption of species coexistence is made, the only L-dependent variables are 




FIG. 5: Daisyworld with additional grey daisy species with albedo 0.5. Other parameters as for 
Fig 4. a) Temperature - thick line shows planetary average, thin lines denote local temperatures 
over daisies. Note that the grey-only solution is stable because it maximises life, despite having 
no temperature regulatory effect, b) Coverage: Thick line denotes bare ground, thin lines denote 
individual species. Note that graph is continuous, with bare ground (x) falling below 0.30625 only 
in the single-species regions. 
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FIG. 6: Bold line shows the mean amount of life from the logistic map, plotted against recursion 
parameter a. Background dots show the bifurcation diagram. A small amount of random noise 
added to y{n) (drawn from a flat distribution ±0.001) eliminates most of the ordered behaviour, 
and the gives mean population above 0.5 throughout. Without noise, the mean population is still 
above 0.5 for all the chaotic cases, and for all the ordered windows except for the tiny region between 
3.9602 and 3.9616, 0.04% of the total. Populations above 0.5 are found for stable attractors across 
the family of logistic maps of the form y„+i = a[j/n(l — VnW 
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